*Code for estimating margin effects from IHS-linear dummy variable coefficients
*Following Bellemare and Wichman equation (10)


** Early output used to create Table 1 in the paper (but table was typed up, rather than output)
** Creates Appendix Table F1

capture log close
log using reg_mfx, replace
timer on 1


************************************
*create program to generate value of xb+ue
program calc_xb_ue

*ue is avg value of the sum of cell fe + idiosyncratic error.
predict ue, ue
summarize ue, meanonly
scalar uebar=r(mean)

*xb needs to sum all continent*year betas times avg value of each continent*year value, plus same for temperature controls. Needs to be a scalar.
*cont*yr effects

*generate means of continent FEs
forvalues c=1/5 {
sum i`c'.continent, meanonly
sca cont`c'bar=r(mean) 
}

* generate means of continentxyear FEs
* for some reason, Stata excludes 2000 for Africa and 2011 for other continents so use Africa 2011 to start summation
* fully balanced panel so average value for each yr dummy is 1/12
gen contyr_atmeans=_b[1.continent#2011.year]*cont1bar*(1/12)
forvalues y=2000/2010 {
	forvalues c=1/5{
	replace contyr_atmeans=contyr_atmeans+ _b[`c'.continent#`y'.year]*cont`c'bar*(1/12)
}	
}
*contyr_atmeans is a same for all obs so just take first value 
sca cy_atmeans=contyr_atmeans[1] 

*cubic in temp dev
sum dev_tmp, meanonly
sca dev_tmpbar=r(mean)
sca temp_atmeans = _b[dev_tmp]*dev_tmpbar + _b[c.dev_tmp#c.dev_tmp]*dev_tmpbar^2 + _b[c.dev_tmp#c.dev_tmp#c.dev_tmp]*dev_tmpbar^3


*then sum mean(ue)+scalar xb. Sum of constant, cont*yr effects, temp effects, cell FE, and error. (Everything but drought var.) This scalar will be what is in the parens after sinh in the denominator.
sca contyr_temp_cons_ue=_b[_cons]+cy_atmeans+temp_atmeans+uebar

*need to drop variables before next use
drop ue contyr_atmeans

end
*end of program, now deploy for various models
*********************************************************

use for_reg.dta, clear

egen cellid=group(x y)
xtset cellid year

**# Table 1 in paper: binary drought coefficients with marginal effects
foreach di in dsi pdsi {
if "`di'"=="dsi" {
local modeln "1"
}
else{
local modeln "2"	
}

xtreg ihslights esm`di' c.dev_tmp##c.dev_tmp##c.dev_tmp i.continent#i.year if flare==0, fe vce(cluster conpfaf4) 

quietly: calc_xb_ue

collect get beta=(_b[esm`di']), tags(model[`modeln'] coeff[Drought])
collect get beta_se=(_se[esm`di']), tags(model[`modeln'] coeff[Drought])
collect, tag(model[`modeln'] coeff[Drought]): nlcom (sinh(_b[esm`di']+contyr_temp_cons_ue)/sinh(contyr_temp_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esm`di']+contyr_temp_cons_ue)/sinh(contyr_temp_cons_ue))-1)), tags(model[`modeln'] coeff[Drought])

sca drop _all
}


collect label levels result gdp "Predicted effect on GDP"
collect label levels model 1 "DSI" 2 "sc-PDSI", modify
collect label levels result _r_b "Marginal effect", modify
collect label levels result beta "Coefficient on moderate or worse drought", modify
collect label levels result beta_se "", modify
collect style cell result[_r_b], nformat(%6.4f) 
collect style cell result[_r_se], sformat((%s)) nformat(%6.4f) 
collect style cell result[beta], nformat(%6.4f) 
collect style cell result[beta_se], sformat((%s)) nformat(%6.4f) 
collect style cell result[gdp], nformat(%6.4f)  
collect stars _r_p 0.01 "**" 0.05 "*" 0.1 "+", attach(_r_b)
collect layout (result[beta beta_se] result[_r_b _r_se] result[gdp]) (model)
collect style putdocx, layout(autofitcontents) title(Binary drought effects) 

putdocx begin
putdocx collect

putdocx save "../output/table1.docx", replace



**# Mfx for Table 2 (GW), column 2 model 


gen alluv=(aqtyp_max=="Major Alluvial")
gen loc_shal=(aqtyp_max=="Local/Shallow")
gen complex=(aqtyp_max=="Complex")
gen karst=(aqtyp_max=="Karst")

label var alluv "Major alluvial"
label var loc_shal "Local/shallow"
label var complex "Complex"
label var karst "Karst"


*interactions with drought
foreach var in alluv loc_shal complex karst  {
  local lbl: variable label `var'
  gen esmdsi_`var'=esmdsi*`var'
   label var esmdsi_`var' "`lbl'*drought"
 }	

xtreg ihslights esmdsi esmdsi_alluv esmdsi_loc_shal esmdsi_complex c.dev_tmp##c.dev_tmp##c.dev_tmp i.continent#i.year if flare==0 & !missing(aqtyp_max), fe vce(cluster conpfaf4)

quietly: calc_xb_ue

*groundwater-drought interactions means
foreach type in alluv loc_shal complex {
sum esmdsi_`type', meanonly
sca esmdsi_`type'bar=r(mean)
}
sum esmdsi, meanonly
sca esmdsibar=r(mean)

*FIRST COEFF (esmdsi)
sca gw_intx = _b[esmdsi_alluv]*esmdsi_alluvbar + _b[esmdsi_loc_shal]*esmdsi_loc_shalbar + _b[esmdsi_complex]*esmdsi_complexbar

sca contyr_temp_gwintx_cons_ue=gw_intx+contyr_temp_cons_ue

di "Table 2, column 2 model"
di "For comparison _b[esmdsi] is " _b[esmdsi]
collect get beta=(_b[esmdsi]), tags(model[4] coeff[Drought])
collect, tag(model[4] coeff[Drought]): nlcom (sinh(_b[esmdsi]+contyr_temp_gwintx_cons_ue)/sinh(contyr_temp_gwintx_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi]+contyr_temp_gwintx_cons_ue)/sinh(contyr_temp_gwintx_cons_ue))-1)), tag(model[4] coeff[Drought])

*SECOND COEFF (esmdsi*alluv)

sca gw_vars = _b[esmdsi]*esmdsibar + _b[esmdsi_loc_shal]*esmdsi_loc_shalbar + _b[esmdsi_complex]*esmdsi_complexbar

sca contyr_temp_gwvars_cons_ue=gw_vars+contyr_temp_cons_ue

di "Table 2, column 2 model"
di "For comparison _b[esmdsi_alluv] is " _b[esmdsi_alluv]
collect get beta=(_b[esmdsi_alluv]), tags(model[4] coeff[Drought*alluv])
collect, tag(model[4] coeff[Drought*alluv]): nlcom (sinh(_b[esmdsi_alluv]+contyr_temp_gwvars_cons_ue)/sinh(contyr_temp_gwvars_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_alluv]+contyr_temp_gwvars_cons_ue)/sinh(contyr_temp_gwvars_cons_ue))-1)), tag(model[4] coeff[Drought*alluv])

sca drop _all

**# Mfx for Table 4, column 1 model 

* characteristics for interactions
gen nonhydro=((count_dams-count_hydro)>0)
replace rescap=rescap/10^3
replace rescaphydro=rescaphydro/10^3
gen up1_nonhydro=((up1_count_dams-up1_count_hydro)>0)
gen farupdam=((upstream_dams-up1_count_dams)>0)
gen pop=pop2000/10^3

label var nonhydro "Non-hydro dam"
label var farupdam "Far upstream dam"
label var pop "Population density"
label var up1_nonhydro "Near upstream non-hydro"

*explicit interactions for cleaner output
foreach var in ifdam ifhydro nonhydro pop cti rescap rescaphydro up1_ifdam up1_ifhydro up1_nonhydro farupdam {
	local lbl: variable label `var'
*extreme/severe/moderate
	gen esmdsi_`var'=(dsi_cat11<=3)*`var'
	label var esmdsi_`var' "`lbl'*drought"
	}


quietly: xtreg ihslights esmdsi esmdsi_ifdam c.dev_tmp##c.dev_tmp##c.dev_tmp i.year#i.continent if flare==0, fe vce(cluster conpfaf4)

quietly: calc_xb_ue

*dam-drought interaction means
sum esmdsi_ifdam, meanonly
sca esmdsi_ifdambar=r(mean)
sum esmdsi, meanonly
sca esmdsibar=r(mean)


*FIRST COEFF (esmdsi)
*then sum mean(ue)+scalar xb. Sum of constant, cont*yr effects, temp effects, cell FE, and error. (Everything but drought var.) This scalar will be what is in the parens after sinh in the denominator.

sca dam_intx = _b[esmdsi_ifdam]*esmdsi_ifdambar
sca contyr_temp_damintx_cons_ue=dam_intx+contyr_temp_cons_ue

*then do nlcom (sinh(scalarxb+ue as number + b_esmdsi)/sinh(scalar xb+ue as number)) - 1

di "*Table 4, column 1 model "
di "for comparison _b[esmdsi] is" _b[esmdsi]
collect get beta=(_b[esmdsi]), tags(model[5] coeff[Drought])
collect, tag(model[5] coeff[Drought]): nlcom (sinh(_b[esmdsi]+contyr_temp_damintx_cons_ue)/sinh(contyr_temp_damintx_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi]+contyr_temp_damintx_cons_ue)/sinh(contyr_temp_damintx_cons_ue))-1)), tag(model[5] coeff[Drought])

*SECOND COEFF (esmdsi*ifdam)
sca drought_var = _b[esmdsi]*esmdsibar
sca contyr_temp_droughtvar_cons_ue=drought_var+contyr_temp_cons_ue

di "*Table 4, column 1 model: dam interaction "
di "for comparison _b[esmdsi_ifdam] is" _b[esmdsi_ifdam]
collect get beta=(_b[esmdsi_ifdam]), tags(model[5] coeff[Drought*dam])
collect, tag(model[5] coeff[Drought*dam]): nlcom (sinh(_b[esmdsi_ifdam]+contyr_temp_droughtvar_cons_ue)/sinh(contyr_temp_droughtvar_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifdam]+contyr_temp_droughtvar_cons_ue)/sinh(contyr_temp_droughtvar_cons_ue))-1)), tag(model[5] coeff[Drought*dam])

sca drop _all


**# Mfx for Table 4, column 2 model (hydro interaction)

quietly: xtreg ihslights esmdsi esmdsi_ifdam esmdsi_ifhydro c.dev_tmp##c.dev_tmp##c.dev_tmp i.year#i.continent if flare==0, fe vce(cluster conpfaf4)

quietly: calc_xb_ue

*dam-drought interactions
foreach v in ifdam ifhydro {
sum esmdsi_`v', meanonly
sca esmdsi_`v'bar=r(mean)
}
sum esmdsi, meanonly
sca esmdsibar=r(mean)

*FIRST COEFF (esmdsi)
sca dam_intx = _b[esmdsi_ifdam]*esmdsi_ifdambar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar

sca contyr_temp_damintx_cons_ue=dam_intx+contyr_temp_cons_ue


di "*Table 4, column 2 model"
di "for comparison _b[esmdsi] is" _b[esmdsi]
collect get beta=(_b[esmdsi]), tags(model[6] coeff[Drought])
collect, tag(model[6] coeff[Drought]): nlcom (sinh(_b[esmdsi]+contyr_temp_damintx_cons_ue)/sinh(contyr_temp_damintx_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi]+contyr_temp_damintx_cons_ue)/sinh(contyr_temp_damintx_cons_ue))-1)), tag(model[6] coeff[Drought])

*SECOND COEFF (esmdsi*ifdam)

sca drought_hydro = _b[esmdsi]*esmdsibar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar
sca contyr_temp_droughthydro_cons_ue=drought_hydro+contyr_temp_cons_ue


di "*Table 4, column 2 model: dam interaction "
di "for comparison _b[esmdsi_ifdam] is" _b[esmdsi_ifdam]
collect get beta=(_b[esmdsi_ifdam]), tags(model[6] coeff[Drought*dam])
collect, tag(model[6] coeff[Drought*dam]): nlcom (sinh(_b[esmdsi_ifdam]+contyr_temp_droughthydro_cons_ue)/sinh(contyr_temp_droughthydro_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifdam]+contyr_temp_droughthydro_cons_ue)/sinh(contyr_temp_droughthydro_cons_ue))-1)), tag(model[6] coeff[Drought*dam])


*THIRD COEFF (esmdsi*ifhydro)

sca drought_dam = _b[esmdsi]*esmdsibar + _b[esmdsi_ifdam]*esmdsi_ifdambar
sca contyr_temp_droughtdam_cons_ue=drought_dam+contyr_temp_cons_ue

di "*Table 4, column 2 model: hydro interaction "
di "for comparison _b[esmdsi_ifhydro] is" _b[esmdsi_ifhydro]
collect get beta=(_b[esmdsi_ifhydro]), tags(model[6] coeff[Hydro*drought])
collect, tag(model[6] coeff[Hydro*drought]):  nlcom (sinh(_b[esmdsi_ifhydro]+contyr_temp_droughtdam_cons_ue)/sinh(contyr_temp_droughtdam_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifhydro]+contyr_temp_droughtdam_cons_ue)/sinh(contyr_temp_droughtdam_cons_ue))-1)), tag(model[6] coeff[Hydro*drought])

**# Mfx for Table 4, column 3 model 
xtreg ihslights esmdsi esmdsi_ifdam esmdsi_ifhydro esmdsi_pop c.dev_tmp##c.dev_tmp##c.dev_tmp i.year#i.continent if flare==0, fe vce(cluster conpfaf4)

quietly: calc_xb_ue

*dam-drought interactions
foreach v in ifdam ifhydro pop {
sum esmdsi_`v', meanonly
sca esmdsi_`v'bar=r(mean)
}
sum esmdsi, meanonly
sca esmdsibar=r(mean)

*FIRST COEFF (esmdsi)

sca dam_intx_pop = _b[esmdsi_ifdam]*esmdsi_ifdambar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar + _b[esmdsi_pop]*esmdsi_popbar

sca contyr_temp_damintpop_cons_ue=dam_intx_pop+contyr_temp_cons_ue

di _b[esmdsi]
collect get beta=(_b[esmdsi]), tags(model[7] coeff[Drought])
collect, tag(model[7] coeff[Drought]): nlcom (sinh(_b[esmdsi]+contyr_temp_damintpop_cons_ue)/sinh(contyr_temp_damintpop_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi]+contyr_temp_damintpop_cons_ue)/sinh(contyr_temp_damintpop_cons_ue))-1)), tag(model[7] coeff[Drought])

*SECOND COEFF (esmdsi*ifdam)

sca drought_hydro_pop = _b[esmdsi]*esmdsibar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar + _b[esmdsi_pop]*esmdsi_popbar
sca contyr_temp_drthydpop_cons_ue=drought_hydro_pop+contyr_temp_cons_ue

di "_b[esmdsi_ifdam] is " _b[esmdsi_ifdam]
collect get beta=(_b[esmdsi_ifdam]), tags(model[7] coeff[Drought*dam])
collect, tag(model[7] coeff[Drought*dam]): nlcom (sinh(_b[esmdsi_ifdam]+contyr_temp_drthydpop_cons_ue)/sinh(contyr_temp_drthydpop_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifdam]+contyr_temp_drthydpop_cons_ue)/sinh(contyr_temp_drthydpop_cons_ue))-1)), tag(model[7] coeff[Drought*dam])

*THIRD COEFF (esmdsi*ifhydro)

sca drought_dam_pop = _b[esmdsi]*esmdsibar + _b[esmdsi_ifdam]*esmdsi_ifdambar + _b[esmdsi_pop]*esmdsi_popbar
sca contyr_temp_drtdampop_cons_ue=drought_dam_pop+contyr_temp_cons_ue

di "_b[esmdsi_ifhydro] is " _b[esmdsi_ifhydro]
collect get beta=(_b[esmdsi_ifhydro]), tags(model[7] coeff[Hydro*drought])
collect, tag(model[7] coeff[Hydro*drought]): nlcom (sinh(_b[esmdsi_ifhydro]+contyr_temp_drtdampop_cons_ue)/sinh(contyr_temp_drtdampop_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifhydro]+contyr_temp_drtdampop_cons_ue)/sinh(contyr_temp_drtdampop_cons_ue))-1)), tag(model[7] coeff[Hydro*drought])

*Table 4, column 4 model 

xtreg ihslights esmdsi esmdsi_ifdam esmdsi_rescap esmdsi_ifhydro esmdsi_rescap esmdsi_rescaphydro c.dev_tmp##c.dev_tmp##c.dev_tmp i.year#i.continent if flare==0, fe vce(cluster conpfaf4)

*dam-drought interaction means
foreach v in ifdam ifhydro rescap rescaphydro  {
sum esmdsi_`v', meanonly
sca esmdsi_`v'bar=r(mean)
}
sum esmdsi, meanonly
sca esmdsibar=r(mean)

*FIRST COEFF (esmdsi)
sca dam_intx_rescap = _b[esmdsi_ifdam]*esmdsi_ifdambar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar + _b[esmdsi_rescap]*esmdsi_rescapbar + _b[esmdsi_rescaphydro]*esmdsi_rescaphydrobar

sca contyr_temp_damintres_cons_ue=dam_intx_rescap+contyr_temp_cons_ue

di _b[esmdsi]
collect get beta=(_b[esmdsi]), tags(model[8] coeff[Drought])
collect, tag(model[8] coeff[Drought]): nlcom (sinh(_b[esmdsi]+contyr_temp_damintres_cons_ue)/sinh(contyr_temp_damintres_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi]+contyr_temp_damintres_cons_ue)/sinh(contyr_temp_damintres_cons_ue))-1)), tag(model[8] coeff[Drought])

*SECOND COEFF (esmdsi*ifdam)
sca drought_hydro_rescap = _b[esmdsi]*esmdsibar + _b[esmdsi_ifhydro]*esmdsi_ifhydrobar + _b[esmdsi_rescap]*esmdsi_rescapbar + _b[esmdsi_rescaphydro]*esmdsi_rescaphydrobar
sca contyr_temp_drthydpres_cons_ue=drought_hydro_rescap+contyr_temp_cons_ue
collect get beta=(_b[esmdsi_ifdam]), tags(model[8] coeff[Drought*dam])
collect, tag(model[8] coeff[Drought*dam]): nlcom (sinh(_b[esmdsi_ifdam]+contyr_temp_drthydpres_cons_ue)/sinh(contyr_temp_drthydpres_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifdam]+contyr_temp_drthydpres_cons_ue)/sinh(contyr_temp_drthydpres_cons_ue))-1)), tag(model[8] coeff[Drought*dam])


*THIRD COEFF (esmdsi*ifhydro)

sca drought_dam_res = _b[esmdsi]*esmdsibar + _b[esmdsi_ifdam]*esmdsi_ifdambar + _b[esmdsi_rescap]*esmdsi_rescapbar + _b[esmdsi_rescaphydro]*esmdsi_rescaphydrobar
sca contyr_temp_drtdamres_cons_ue=drought_dam_res+contyr_temp_cons_ue

collect get beta=(_b[esmdsi_ifhydro]), tags(model[8] coeff[Hydro*drought])
collect, tag(model[8] coeff[Hydro*drought]): nlcom (sinh(_b[esmdsi_ifhydro]+contyr_temp_drtdamres_cons_ue)/sinh(contyr_temp_drtdamres_cons_ue))-1
collect get gdp=(.3* ((sinh(_b[esmdsi_ifhydro]+contyr_temp_drtdamres_cons_ue)/sinh(contyr_temp_drtdamres_cons_ue))-1)), tag(model[8] coeff[Hydro*drought])

**# Appendix Table F1 layout
 

collect label levels model 1 "Table 1 (DSI)" 2 "Table 1 (sc_PDSI)" 4 "Table 2, col. 2 (GW)" 5 "Table 4 (dams), col 1" 6 "Table 4, col 2" 7 "Table 4 col 3" 8 "Table 4 col 4", modify
collect label levels result _r_b "Marginal effect", modify
collect label levels result beta "Coefficient", modify
collect style cell result[_r_b], nformat(%6.4f) halign(center)
collect style cell result[_r_se], sformat((%s)) nformat(%6.4f)
collect style cell result[beta], nformat(%6.4f) halign(center)
collect style cell result[gdp], nformat(%6.4f) halign(center)
collect stars _r_p 0.01 "**" 0.05 "*" 0.1 "+", attach(_r_b)
collect layout (model#coeff) (result[beta] result[_r_b _r_se] result[gdp])
collect style putdocx, layout(autofitcontents) title(Marginal effects estimates at mean) 

putdocx begin
putdocx collect

putdocx save "../output/tableF1_mfx.docx", replace




collect clear
putdocx clear

timer off 1
timer list 1
log close

timer clear 1 